import typing
import firedrake
import numpy as np
import ufl
/home/firedrake/firedrake/lib/python3.10/site-packages/pytools/__init__.py:2449: UserWarning: unable to find git revision
warn("unable to find git revision")
firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance
import viskex
Generate and plot meshes
interval = firedrake.UnitIntervalMesh(6)
square_tria = firedrake.UnitSquareMesh(3, 3, quadrilateral=False, diagonal="right")
square_quad = firedrake.UnitSquareMesh(3, 3, quadrilateral=True, diagonal="right")
cube_tetra = firedrake.UnitCubeMesh(3, 3, 3, hexahedral=False)
cube_hexa = firedrake.UnitCubeMesh(3, 3, 3, hexahedral=True)
for mesh in (interval, square_tria, square_quad, cube_tetra, cube_hexa):
mesh.topology_dm.removeLabel(firedrake.cython.dmcommon.FACE_SETS_LABEL)
del mesh
viskex.firedrake.plot_mesh(interval)
viskex.firedrake.plot_mesh(interval, dim=1)
viskex.firedrake.plot_mesh(interval, dim=0)
viskex.firedrake.plot_mesh(square_tria)
viskex.firedrake.plot_mesh(square_tria, dim=2)
viskex.firedrake.plot_mesh(square_tria, dim=1)
viskex.firedrake.plot_mesh(square_tria, dim=0)
viskex.firedrake.plot_mesh(square_quad)
viskex.firedrake.plot_mesh(square_quad, dim=2)
viskex.firedrake.plot_mesh(square_quad, dim=1)
viskex.firedrake.plot_mesh(square_quad, dim=0)
viskex.firedrake.plot_mesh(cube_tetra)
viskex.firedrake.plot_mesh(cube_tetra, dim=3)
viskex.firedrake.plot_mesh(cube_tetra, dim=2)
viskex.firedrake.plot_mesh(cube_tetra, dim=1)
viskex.firedrake.plot_mesh(cube_tetra, dim=0)
viskex.firedrake.plot_mesh(cube_hexa)
viskex.firedrake.plot_mesh(cube_hexa, dim=3)
viskex.firedrake.plot_mesh(cube_hexa, dim=2)
viskex.firedrake.plot_mesh(cube_hexa, dim=1)
viskex.firedrake.plot_mesh(cube_hexa, dim=0)
Generate and plot subdomains
def mark_subdomains(mesh: firedrake.MeshGeometry) -> firedrake.MeshGeometry: # type: ignore[no-any-unimported]
"""Mark left and right subdomains in a given mesh with values 1 and 2, respectively."""
cellname = mesh.ufl_cell().cellname()
if cellname in ("interval", "triangle", "tetrahedron"):
subdomains_function_space = firedrake.FunctionSpace(mesh, "DP", 0)
elif cellname in ("quadrilateral", "hexahedron"):
subdomains_function_space = firedrake.FunctionSpace(mesh, "DQ", 0)
else:
raise RuntimeError("Invalid cellname")
x = firedrake.SpatialCoordinate(mesh)
left_subdomain = firedrake.Function(subdomains_function_space).interpolate(
firedrake.conditional(x[0] <= 1.0 / 3.0, 1.0, 0.0))
right_subdomain = firedrake.Function(subdomains_function_space).interpolate(
firedrake.conditional(x[0] >= 2.0 / 3.0, 1.0, 0.0))
mesh_with_subdomains = firedrake.RelabeledMesh(mesh, [left_subdomain, right_subdomain], [1, 2])
mesh_with_subdomains.init()
return mesh_with_subdomains
interval_with_subdomains = mark_subdomains(interval)
square_tria_with_subdomains = mark_subdomains(square_tria)
square_quad_with_subdomains = mark_subdomains(square_quad)
cube_tetra_with_subdomains = mark_subdomains(cube_tetra)
cube_hexa_with_subdomains = mark_subdomains(cube_hexa)
viskex.firedrake.plot_mesh_entities(
interval_with_subdomains, 1, "subdomains", interval_with_subdomains.cell_subset(2).indices
)
viskex.firedrake.plot_mesh_entities(
interval_with_subdomains, 1, "subdomains", interval_with_subdomains.cell_subset(2).indices,
2 * np.ones_like(interval_with_subdomains.cell_subset(2).indices)
)
viskex.firedrake.plot_mesh_sets(interval_with_subdomains, 1, "subdomains")
viskex.firedrake.plot_mesh_entities(
square_tria_with_subdomains, 2, "subdomains", square_tria_with_subdomains.cell_subset(2).indices
)
viskex.firedrake.plot_mesh_entities(
square_tria_with_subdomains, 2, "subdomains", square_tria_with_subdomains.cell_subset(2).indices,
2 * np.ones_like(square_tria_with_subdomains.cell_subset(2).indices)
)
viskex.firedrake.plot_mesh_sets(square_tria_with_subdomains, 2, "subdomains")
viskex.firedrake.plot_mesh_sets(square_quad_with_subdomains, 2, "subdomains")
viskex.firedrake.plot_mesh_sets(cube_tetra_with_subdomains, 3, "subdomains")
viskex.firedrake.plot_mesh_sets(cube_hexa_with_subdomains, 3, "subdomains")
Generate and plot boundaries
def mark_boundaries(mesh: firedrake.MeshGeometry) -> firedrake.MeshGeometry: # type: ignore[no-any-unimported]
"""
Mark internal and boundary facets in a given mesh with four different values.
Internal facets of left and right subdomains are associated with values 1 and 2, respectively.
Furthermore, boundary facets on the left and right boundaries are associated with values 3 and 4,
respectively.
"""
cellname = mesh.ufl_cell().cellname()
if cellname in ("interval", ):
boundaries_function_space = firedrake.FunctionSpace(mesh, "P", 1)
elif cellname in ("triangle", "quadrilateral", "tetrahedron", "hexahedron"):
boundaries_function_space = firedrake.FunctionSpace(mesh, "HDiv Trace", 0)
else:
raise RuntimeError("Invalid cellname")
x = firedrake.SpatialCoordinate(mesh)
left_boundary = firedrake.Function(boundaries_function_space).interpolate(
firedrake.conditional(abs(x[0] - 0.) < np.sqrt(np.finfo(float).eps), 1.0, 0.0))
right_boundary = firedrake.Function(boundaries_function_space).interpolate(
firedrake.conditional(abs(x[0] - 1.) < np.sqrt(np.finfo(float).eps), 1.0, 0.0))
left_subdomain = firedrake.Function(boundaries_function_space).interpolate(
firedrake.conditional(
firedrake.And(x[0] <= 1.0 / 3.0, abs(x[0] - 0.) > np.sqrt(np.finfo(float).eps)), 1.0, 0.0))
right_subdomain = firedrake.Function(boundaries_function_space).interpolate(
firedrake.conditional(
firedrake.And(x[0] >= 2.0 / 3.0, abs(x[0] - 1.) > np.sqrt(np.finfo(float).eps)), 1.0, 0.0))
mesh_with_boundaries = firedrake.RelabeledMesh(
mesh, [left_boundary, right_boundary, left_subdomain, right_subdomain], [3, 4, 1, 2])
mesh_with_boundaries.init()
return mesh_with_boundaries
interval_with_boundaries = mark_boundaries(interval)
square_tria_with_boundaries = mark_boundaries(square_tria)
square_quad_with_boundaries = mark_boundaries(square_quad)
cube_tetra_with_boundaries = mark_boundaries(cube_tetra)
# cube_hexa_with_boundaries = mark_boundaries(cube_hexa) # HDiv not implemented
viskex.firedrake.plot_mesh_entities(
interval_with_boundaries, 0, "boundaries",
interval_with_boundaries.exterior_facets.measure_set("exterior_facet", 4).indices
)
viskex.firedrake.plot_mesh_entities(
interval_with_boundaries, 0, "boundaries",
interval_with_boundaries.exterior_facets.measure_set("exterior_facet", "everywhere").size
+ interval_with_boundaries.interior_facets.measure_set("interior_facet", 2).indices
)
viskex.firedrake.plot_mesh_sets(interval_with_boundaries, 0, "boundaries")
viskex.firedrake.plot_mesh_entities(
square_tria_with_boundaries, 1, "boundaries",
square_tria_with_boundaries.exterior_facets.measure_set("exterior_facet", 4).indices
)
viskex.firedrake.plot_mesh_entities(
square_tria_with_boundaries, 1, "boundaries",
square_tria_with_boundaries.exterior_facets.measure_set("exterior_facet", "everywhere").size
+ square_tria_with_boundaries.interior_facets.measure_set("interior_facet", 2).indices
)
viskex.firedrake.plot_mesh_sets(square_tria_with_boundaries, 1, "boundaries")
viskex.firedrake.plot_mesh_sets(square_quad_with_boundaries, 1, "boundaries")
viskex.firedrake.plot_mesh_sets(cube_tetra_with_boundaries, 2, "boundaries")
Interpolate and plot scalar functions
def prepare_scalar_field_cases( # type: ignore[no-any-unimported]
mesh: firedrake.Mesh,
expression: typing.Callable[[ufl.core.expr.Expr], ufl.core.expr.Expr]
) -> typing.Tuple[firedrake.Function, typing.Tuple[ufl.core.expr.Expr, ufl.FunctionSpace]]:
"""Prepare scalar field cases."""
scalar_function_space = firedrake.FunctionSpace(mesh, "CG", 2)
scalar_field_ufl = expression(ufl.SpatialCoordinate(mesh))
scalar_field = firedrake.interpolate(scalar_field_ufl, scalar_function_space)
return scalar_field, (scalar_field_ufl, scalar_function_space)
interval_scalar_field, interval_scalar_field_ufl = prepare_scalar_field_cases(
interval, lambda x: x[0]**3)
square_tria_scalar_field, square_tria_scalar_field_ufl = prepare_scalar_field_cases(
square_tria, lambda x: x[0]**3 + x[1]**2)
square_quad_scalar_field, square_quad_scalar_field_ufl = prepare_scalar_field_cases(
square_quad, lambda x: x[0]**3 + x[1]**2)
cube_tetra_scalar_field, cube_tetra_scalar_field_ufl = prepare_scalar_field_cases(
cube_tetra, lambda x: x[0]**3 + x[1]**2 + x[2]**4)
cube_hexa_scalar_field, cube_hexa_scalar_field_ufl = prepare_scalar_field_cases(
cube_hexa, lambda x: x[0]**3 + x[1]**2 + x[2]**4)
viskex.firedrake.plot_scalar_field(interval_scalar_field, "scalar")
viskex.firedrake.plot_scalar_field(interval_scalar_field_ufl, "scalar")
viskex.firedrake.plot_scalar_field(square_tria_scalar_field, "scalar")
viskex.firedrake.plot_scalar_field(square_tria_scalar_field, "scalar", warp_factor=1.0)
viskex.firedrake.plot_scalar_field(square_tria_scalar_field_ufl, "scalar")
viskex.firedrake.plot_scalar_field(square_quad_scalar_field, "scalar")
viskex.firedrake.plot_scalar_field(square_quad_scalar_field, "scalar", warp_factor=1.0)
viskex.firedrake.plot_scalar_field(square_quad_scalar_field_ufl, "scalar")
viskex.firedrake.plot_scalar_field(cube_tetra_scalar_field, "scalar")
viskex.firedrake.plot_scalar_field(cube_tetra_scalar_field_ufl, "scalar")
viskex.firedrake.plot_scalar_field(cube_hexa_scalar_field, "scalar")
viskex.firedrake.plot_scalar_field(cube_hexa_scalar_field_ufl, "scalar")
Interpolate and plot vector functions
def prepare_vector_field_cases( # type: ignore[no-any-unimported]
mesh: firedrake.Mesh,
expression: typing.Callable[[ufl.core.expr.Expr], ufl.core.expr.Expr]
) -> typing.Tuple[firedrake.Function, typing.Tuple[ufl.core.expr.Expr, ufl.FunctionSpace]]:
"""Prepare vector field cases."""
vector_function_space = firedrake.VectorFunctionSpace(mesh, "CG", 2)
vector_field_ufl = ufl.as_vector(expression(ufl.SpatialCoordinate(mesh)))
vector_field = firedrake.interpolate(vector_field_ufl, vector_function_space)
return vector_field, (vector_field_ufl, vector_function_space)
square_tria_vector_field, square_tria_vector_field_ufl = prepare_vector_field_cases(
square_tria, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
square_quad_vector_field, square_quad_vector_field_ufl = prepare_vector_field_cases(
square_quad, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
cube_tetra_vector_field, cube_tetra_vector_field_ufl = prepare_vector_field_cases(
cube_tetra, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))
cube_hexa_vector_field, cube_hexa_vector_field_ufl = prepare_vector_field_cases(
cube_hexa, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))
viskex.firedrake.plot_vector_field(square_tria_vector_field, "vector")
viskex.firedrake.plot_vector_field(square_tria_vector_field, "vector", glyph_factor=0.1)
viskex.firedrake.plot_vector_field(square_tria_vector_field, "vector", warp_factor=1.0)
viskex.firedrake.plot_vector_field(square_tria_vector_field_ufl, "vector")
viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector")
viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector", glyph_factor=0.1)
viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector", glyph_factor=0.1)
viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector", warp_factor=1.0)
viskex.firedrake.plot_vector_field(cube_tetra_vector_field, "vector")
viskex.firedrake.plot_vector_field(cube_tetra_vector_field, "vector", glyph_factor=0.1)
viskex.firedrake.plot_vector_field(cube_tetra_vector_field, "vector", warp_factor=1.0)
viskex.firedrake.plot_vector_field(cube_hexa_vector_field, "vector")
viskex.firedrake.plot_vector_field(cube_hexa_vector_field, "vector", glyph_factor=0.1)
viskex.firedrake.plot_vector_field(cube_hexa_vector_field, "vector", warp_factor=1.0)